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Abstract 

The dynamics of inertial particles in turbulence is modelled and investigated by means 
of direct numerical simulation of an axisymmetrically expanding homogeneous turbulent 
strained flow. This flow can mimic the dynamics of particles close to stagnation points. 
The influence of mean straining flow is explored by varying the dimensionless strain rate 
parameter Sko/eo from 0.2 to 20. We report results relative to the acceleration variances 
and probability density functions for both passive and inertial particles. A high mean strain 
is found to have a significant effect on the acceleration variance both directly, through an 
increase in wave number magnitude, and indirectly, through the coupling of the fluctuat¬ 
ing velocity and the mean flow field. The influence of the strain on normalized particle 
acceleration pdfs is more subtle. For the case of passive particle we can approximate the 
acceleration variance with the aid of rapid distortion theory and obtain good agreement 
with simulation data. For the case of inertial particles we can write a formal expressions 
for the accelerations. The magnitude changes in the inertial particle acceleration variance 
and the effect on the probability density function are then discussed in a wider context for 
comparable flows, where the effects of the mean flow geometry and of the anisotropy at 
the small scales are present. 


1 Introduction 


The motions of small, passive and inertial particles in turbulence has been extensively studied 
in recent years, both from the experimental and theoretical viewpoints. This was motivated by 
a broad range of applications such as the spread of pollutants in the atmosphere and oceans, 
the process of rain and ice f ormati o n in c loud, the transport of sediments i n rivers and estu¬ 
aries, see e.g. the reviews of Shaw ( 2003 1 and Toschi &: Bodenschat d ( 2009ll . Progress in the 
understanding was on the one hand made possible by recent improvements in Lagrangian mea¬ 
surements through particle tracking metho d ologies, result i ng in part from rapid a d vances in 
high speed imaging ( Virant fc Dracos ([i9^;[Ott^]Man3 (l2000ll : IVoth et al\ (l2002h : IXu 
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( 2008h'l. and on t he other hand by inc reased computational capabilities of numerical simulations 


( Yeung fc Pooel (1991); Celani (2007)). 

The objective of this work is to investigate the effects of flow straining on the Lagrangian 
dynamics of small, sub-Kolmogorov scale, passive and inertial particles. Our motivation stems 
both from the fact that many practical turbulent flows are subject to straining motions, such 
as the external flows over bluff or streamlined bodies and internal flows in variable cross sec¬ 


tions (Batchelor b Hunt fc CarruthersI ( 19901: Warhaft ( 198f)l ): Chen et at. 

( 200(ih T lAvvalasomavaiula fc Warhaft (12006 1: Cualtieri fc Meneveaul (|2oiofl l. A mean straining 
flow naturally appears in the proximity of stagnation points. Flow straining is furthermore of 
fundamental interest since it induces a scale dependent anisotropy; the sma llest scales of the flow 
may b e nearly isotropic, whereas the largest scales are highly anisotropic (IBiferale fc Procaccia 
(1200, TM . 

Furthermore, many flows naturally combine straining geometries and inertial particles. The 
flow geometry presented here, namely particle laden turbulent flow undergoing a n axisymmetric 
expansion, has similarities with combustor diffusers in jet engines (jKleinl (Il995r) l. where liquid 
fuel is injected in an expanding flow, and with the flow in the combustion chamber in an intern al 
combustion engine during the compression stroke of the fuel air mixture ( Han fc Reit d ( 199, 5ll ). 

While significant attention has been given to the study of Lagrangian acceleration statistics 
in isotropic turbulence, less attention has been paid to the implications of anisotropic large 
scale flow geometry on the Lagrangian dynamics. Recent experimental and numerical work on 
the Lagrangian behavior of inertial particles in shear flows and i n turbulent bounda r y laye r 


has shown pr o nounc e d effects on th e inertial pa rticle statistics (jGerashchenco et al\ (l2008fl : 


Lavezzo et all ( 2011)1) : Gualtieri et all ( 20091 [201^ ). The persistent small scale anisotropy has 
been found to influence the geometry and alignment of particle clusters and relative particle pair 
velocities. In addition, the combined effects of gravity and shear on particle acceleration variance 
results in an increase in magnitude with the Stokes number. As a consequence, the acceleration 
probability distribution functions (pdfs) became increasingly narrow and skewed with inertia. 
Here, we address a related topic, namely the complexity introduced in the Lagrangian dynamics 
of tracer and inertial particles due to flow straining. 

In an effort to realize the effects of anisotropy in the particle dynamics, we numerically 
simulate axisymmetric expansion of initially isotropic turbulence. The flow is seeded with 
infinitesimal tracer and inertial particles of varied Stokes numbers. We measure the particle 
velocity and acceleration statistics, including variances and probability density functions for 
different strain rates and Stokes numbers. Comparisons are made with predictions of rapid 
distortion theory on tracer accelerations, and the solutions of the Stokes equations for inertial 
particles in the straining flow. 

The paper is organized as follows. In section 2 we briefly introduce the numerical methods for 
simulating an axisymmetric turbulence and particle movements. Parameters of the simulations 
are listed. Section 3 presents the underlying flow field. We discuss our main findings of particle 
acceleration variances and probability density functions in simulation data with the support of 
theoretical estimations in section 4. We also discuss our result in the context of previous work 
in shear flows. In section 5 we present our conclusions. 
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Figure 1: A sketch of the deformation of the simulation domain under straining. The mean flow 
U = {-2Sx, Sy, Sz) corresponds to an ideal flow onto a flat plate. The deforming domain is 
initially elongated in the x-direction but becomes wider in y and z directions with time. Arrows 
indicate the directions of the stream lines of the induced mean flow. 


2 Methodology 

2.1 Flow equations and flow simulation 

The equations describing the motion of an incompressible Newtonian fluid are the continuity 
equation and the Navier-Stokes equations, respectively: 


= 0 , 

= vV'^u. ( 1 ) 

Here, u and p are the instantaneous flow velocity and pressure, respectively, and v is the 
kinematic viscosity of the fluid. 

In this paper we are concerned with a turbulent fluid undergoing an axisymmetric expansion, 
where the mean flow field is described by 

U={-2Sx,Sy,Sz), (2) 

here S is the mean strain rate S = ^ ^ , and Sij = ^ is the mean rate 

of strain tensor. The mean flow corresponds to an ideal flow onto a flat plate, and is realized 
in the flow between contracting pistons or in the expanding flow through a diffuser. Figure [T] 
shows a sketch of the mean flow field, namely streamlines of the mean field, and the deforming 
domain. 

Applying the Reynolds decomposition, and expressing ([1]) in terms of the vector potential 
6, with M = V X 6, one obtains 

- dtV'^b - V X (m X w) + 2Sxd^y^b - SydyV^b - Szd^V^b + SV^b - iSV^biii = (3) 

where u = ii-U is the velocity fluctuation, ui is the vorticity, defined as uj = V x u = -V^b. 
bi is the first component of b and ei is the unit vector in the x-direction. In the following we 
briefly outline the numerical algorithm. 


V • u 
du 

+ u ■ S/u + S/p 
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As in ICvlfason et all (12011 ) in o rder to solve ([3|) we apply a pseudo-spectral method with 
Rogallo’s algorithm ( Rogallol (1981)) where the following variable transformations are per¬ 
formed: 

t' = t. 


X = 


y' = e 


z' = 


and hence the vector potential of the velocity fluctuation satisfies 


dt'S/ "^b - S/' X {u X uj) + SV^b - 3S'v'^5iei = 


-W'^b, 


( 4 ) 


where V' = [e'^^*dx',e~^*dy',e~^*dz'). By adopting this new coordinate system, the physical 
domain deforms with time while the computational lattice grid is time independent, and the 
flow equations become periodic. We then apply the pseudo-s pectral method to equ ations in 
More details about the numerical method can be found in lGvlfason et all (|201l[ ). 


Numerical simulations of this axisymmetric expansion flow were carried out on a Cartesian 
grid with 1024x256x256 and 2048x512x512 grid points in the x, y and z directions, respectively. 
The initial configurations are derived from statistically independent homogeneous and isotropic 
flow simulations which have reached a stationary state after more than 5 large-eddy turnover 
times. The Reynolds numbers, based on the Taylor microscale Aq were i?Ao = 117 and R\q = 193 
before the straining was applied, for the lower and higher grid resolution respectively. Initially 
the physical domain size was [0,87r] x [0,27r] x [0,27r] in x, y and 2 directions respectively, and 
the simulation was terminated when the domain has reached [0,1. Itt] x [0,5.47r] x [0,5.47r] to 
prevent the physical domain to become too flattened. 

The top of Figure [2] shows snapshots of the fluctuating velocity magnitude at three time 
instants during the straining. From left to right, the non-dimensional times are S x t = 0.08 
(shortly after the mean strain is applied), a.t S x t = 0.64, and at S' x t = 0.96 (just before the 
strain simulation is terminated due the large deformation of the physical domain). Additionally, 
the figure shows the coordinate system adopted in the text, and the geometry of the simulation 
domain selected and its deformation. Production of turbulence overwhelms dissipation during 
the straining, reflected in an increase in the turbulent kinetic energy, most notably in the 
compressed component (x). This can be seen in the warmer colors in the rightmost plot. The 
bottom of Figure[2]shows isosurfaces of non-dimensional vorticity w/= 3.17 at the same 
time instants as above, which is respectively 4.36, 2.86 and 2.43 standard deviations above the 
mean vorticity magnitude at the three time instants. From the increased number and size of 
the filaments, we observe that the vorticity is intensified during straining, and the filaments are 
found to gradually align with the y, z-plane due to the mean flow extension in the plane. 


2.2 Equations for particle movements 

To study Lagrangian aspects of this flow we seed the flow with tracers and inertial particles. 
Here, we are concerned with particles that are small compared to the smallest length scales 
present in the flow and their densities considerably higher than the fluid density. The particle 
number densities are furthermore assumed to be sufficiently low so that particle-particle inter¬ 
actions can be ignored. Under the above approximations, the coupling of the particles with the 
carrier fluid can be ignored. 

The Lagrangian equations of inertial particle motion is derived from Newton’s second law, 
and represents the balance between the forces acting on the particles (inertia and Stokes drag). 
The equations describing the motion of a particle of diameter dp and density Pp, located at Xp 
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Figure 2: Snapshots of the magnitude of the fluctuating flow velocity (top) and isosurface of 
the magnitude of vorticity a; = ||V x m|| (bottom) in the deforming domain in a realization of 
the axisymmetrically expanding flow. Left to right: from the onset of the straining simulation 
to the end of the straining simulation at time instants S xt = 0.08,0.64, and 0.96. The size of 
the simulation domain is noted in the figure, and its deformation is displayed. The coordinate 
system adopted in the text is also shown. We simulated the axisymmetric turbulence with two 
resolutions using 2048 x 512 x 512 and 1024 x 256 x 256 computational nodes. The illustration 
shows one of the realizations of the 1024 x 256 x 256 simulations with S = 10. The isosurfaces 
plotted have non-dimensional vorticity = 3.17, which is respectively 4.36, 2.86 and 

2.43 standard deviations above the mean vorticity magnitude at the three time instants and 
chosen to illustrate the flow structure. 
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and with instantaneous velocity Vp are 

dxp 

dt 

dvp 1 . 

— = —{u{Xp)-Vp), 

CLl Tp 

where Tp = [ddpIlSv is the Stokes relaxation time for the particle and (3 = {pp - Pf)/pf is 
the relative density ratio between the particle and the fluid. The Stokes number St = TpjTr) 
characterizes the inertia of a particle in the flow, where is the Kolmogorov timescale of the 
flow. 

For the tracer particles (zero inertia), the particle velocity is the same as the fluid velocity 
at the particle location, which yields the kinematic relation: 

The ordinary differential equations ([5]), ([5]), and © are solved numerically by the second 
order Adams-Bashforth method. In equations (jB]) and the instantaneous flow velocity at 
the particle location, Xp, is evaluated as 

~ ^ (^p) 5 

that is, the mean flow velocity is evaluated at the location of the particle through the formula 
U{xp) = (-2Sxp, Sup, Szp), and the flow velocity fluctuation is interpolated to the particle 
position. 

We initialize the particles and the fluid velocity with steady state homogeneous isotropic 
simulations. The particles are uniformly distributed over the domain prior to the forced homo¬ 
geneous isotropic simulation is carried out. At the beginning of the straining, at t = 0, we add 
the mean flow velocity and acceleration component due to the strain geometry to the existing 
particle velocity. We conduct simulations with 1024x256^ and 2048x512^ collocation points. 
For the lower resolution simulation we use 16 independent flow realizations with 5x10^ particles 
of each type (6 different Stokes numbers) and for the higher resolution simulations we perform 
10 independent flow realizations with 4x10^ particles of each type. 

Table[T]shows the various flow parameters for the simulations performed. The range of strain 
rates selected is such that its effect on the smallest scales of the flow range from being negligible 
to substantial. The higher strain rates are felt intensely by the large scale flow, whereas the 
lower strain rates have mild effect on the large scales. The value of the strain parameters, S'r^o 
and Skojeo, which compare the strain time with the local timescales of the flow, indicate the 
importance of the various terms in the evolution equation of the velocity field. 


(5) 

( 6 ) 


IMaxev fc Rilev (1983); Bee et al. (l2006h'l: 


3 Underlying flow fleld 

The left hand side of Figure [3] shows the evolution of the Reynolds stresses {{uf)) normalized 
by the initial turbulent kinetic energy (fco)- The component along the compressed direction 
(xi) grows rapidly in most cases, whereas the components along the expanding directions are 
suppressed or remain roughly constant. At the lowest strain rate all component are suppressed 
during the simulation time. For large times, rapid distortion theory (RDT) predicts an expo¬ 
nential growth of the Reynolds stresses, in the proportions (wf) = 2(Mf) = 2{u\). The kinetic 
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Table 1: Flow parameters in the Direct Numerical Simulations, based on the homogeneous 
isotropic simulation prior to the application of the straining. Here = | (('^lo) + ('*^ 20 ) ('*^ 30 )) 

is the turbulent kinetic energy, eo = is the energy dissipation rate, £0 = 

is the integral length scale, and the Kolmogorov length scale 770 = The subscript 0 

indicates the parameter values are taken prior to the straining. (Units are arbitrary) 


Simulation domain 
Rxo = 

ko=^ (Ko> + (ulo) + (U30» 
Ao/?7o = jrio 

S 

Tr,o = 

4/770 ^ 

_ 1_ / Oujo 9«io. \ 

""O - 2 ' dxj dxj I 


R'^riO 

Skolep 


1024x256x256 

117 

4.6 

20.9 

0.1,0.5,1,4,10 

0.051 

1.75 

164.4 

0.0163 ± 0.0006 
2.18 ± 0.15 
0.0052 

0,0.2,0.3,0.5,1,2 
0.0051,0.0255,0.051,0.204,0.51 
0.21,1.06,2.11,8.44,21.1 


2048x512x512 

193 

4.9 

27.4 

1,4,10 

0.031 

1.81 

332.6 

0.008 ± 0.0004 
2.12 ± 0.4 
0.00205 

0,0.23,0.34,0.56,1.12,2.25 

0.031,0.124,0.31 

2.31,9.25,23.1 


energy increases with time for all the strain rates, although the lower rates display an initial 
drop in energy and a subsequent long term increase. The right hand side of the figure shows the 
evolution of the anisotropy tensor with time. The curves corresponding to the 

lowest strain rates markedly deviate from the others as the turbulent kinetic energy decreases 
during the straining, and the straining motions are fairly mild, even for the largest scales of 
motions. 

The s hort t e rm R DT prediction plotted in Figure [3] is derived from the Reynolds stress 
equation ( Pope ( 2000ll l 


— {u^Uj) = Vi 


+ 


( 8 ) 


d’’) - 


whereVij = -{uiUk)^^-{ujUk)^^ is the production rate of Reynolds stress, + 

is the rapid pressure-rate-of-strain tensor, and is the rapid pressure that satisfies ^ = 

flight before the straining starts, the initial configuration of flow is isotropic, and 
■f’Py ( Pope ( 200r)ll 'l. In an axisymmetric expansion flow, the production rates are 
Vii = 45'((ui)^),P22 = -25'((m 2)^), and V 33 = - 2 S{{u 3 Y). Therefore, in early times RDT 
predicts 




= ((wio)">et''*, ((772)2) = ((7720)">e-^''*, ((773)2) = ( (773o)2)e-1 , (9) 

where (( 77 ^ 0 ) 2 ) represent the initial values of the Reynolds stresses {{{uiY),i = 1,2,3). However, 
in order for rapid distortion theory to apply, the parameters must satisfy » 1 and Sk/e » 1 
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Figure 3: Left: Normalized Reynolds stresses {u^)lko vs. strain time S x t. Right: The flow 

{UiUj) 1 

anisotropy tensor bij = - -r- dij vs. strain time. Solid symbols represent the i = 1 component 

y'Ui'tli/ 3 

and open symbols represent the i = 2 component, fco is the initial kinetic energy of the fluid 
prior to the straining. The diamond (O), pentagram (☆), circle (O), square (□) and triangle 
(A) symbols represent data from S = 0.1,0.5,1,4,10 in the Rxo = 117 How, respectively. Blue, 
green and red lines, indicate the data from S = 1,4,10 in the R\q = 193 flow; solid and dashed 
lines represent the i= \ and i = 2 components. Black solid and dashed lines represent the short 
term RDT predictions ([HI). Estimates of statistical error bars on data from S = 0.1,1,10 are 
shown and are computed according to (nni) with Xj being {uf)lkQ and bu in the j-ih realization. 
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Figure 4: Ratios of the variances of velocity derivatives {{dxU 2 )^)I{{dxUi)^) vs. normalized time 
t/r^o- The circle (O), square (□) and triangle (A) symbols represent strain rates 5 = 1,4,10 from 
the R\o = 117 set, and crosses (x), asterisks (*) and pluses (+) represent strain rates 5 = 1,4,10 
from the R\o = 193 set. An estimate of statistical error bar is shown (within the symbols) and 
computed according to (fTOll with Xj being {{dxU 2 )^)I{{dxUi)^) in the j-th realization. The 
solid line shows the theoretical prediction for this ratio in the isotropic turbulence. 


I Batchelor (195^). Only at the highest rate of strain is the latter constraint weakly satisfied, 
and therefore one does not observe close matches between the predictions of RDT and the 
Reynolds stresses in simulation data. The global anisotropy is much less sensitive, and short 
term RDT predicts the anisotropy well. 

The error bars in Figure [3] indicate the statistical error of the quantities estimated from the 
finite number of realizations of the flow. That is, in X realizations of the turbulent flow with a 
particular strain rate S, one obtains samples {Ai, A 2 ,..., of a quantity A. The estimated 
standard error is the sample standard deviation of {Ai, A 2 ,..., Avr} divided by \/jV. That is. 




( 10 ) 


for data from A realizations; here A is the mean of {Xj}^^ . In this work, the length of the 
symmetric error bars in Figures [3] to [8] is twice of the estimated statistical error. 

Since tracer and inertial particle accelerations occur primarily at the smallest scales of 
motions, it is useful to look at the effects of the straining on the small scales. Figure S] shows 
a measure of the small scale anisotropy, the ratio of variances of longitudinal derivatives of the 
transverse and longitudinal velocity components ((dxU 2 )^)/((dxUi)^) with respect to time. At 
the highest strain rates, the anisotropy due to the straining is present at the smallest scales of 
motion, whereas for the lower strain rates, the flow appears nearly isotropic at the small scales. 
Note that the isotropic prediction for the ratio is ((dxU 2 )^}/((dxUi)^} = 2. The small scale 
anisotropy appears to become close to a constant after an initial transition period for the lower 
strain rates, but a stationary state is not reached at the highest rate of strain for this quantity. 

The left hand side of Figure [5] shows the time evolution of longitudinal derivative skew- 
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Figure 5: Left: Skewness in the compression and expansion directions. Solid lines (—): 
{{dxUi)^)l{{dxUiY)^l'^. Dashed lines (—): {{dyU 2 )^)j{{dyU 2 Y)^^'^■ Right: Kurtosis in the com¬ 
pression and expansion directions. Solid lines (—): {{d^uiY)j{{dxUiY)^. Dashed lines (—): 
{{dyU 2 Y)!{{dyU 2 Y)^■ In both plots: circles (O), squares (□) and triangles (A): S = 1,4,10 
from the R\q = 117 set, and crosses (x), asterisks (=t=) and pluses (+): S'=1,4,10 from the 
i?Ao = 193 set. An estimate of statistical error bar is shown and computed according to (ITTIll 
with Xj being the skewness and kurtosis in the j-th realization. 


ness, along the directions of compression and expansion in the flow. Before the straining is 
applied, the skewness has a value of about [-0.4,-0.5] as expected, but the straining causes 
a marked change in its value and becomes positive in the expanding direction. The effect in 
the compressed direction is more subtle, but an increase in magnitude (larger negative values) 
appear to occur for all strain rates given that the simulation is run for a sufficiently long time. 
The sign change indicates a change in the small scale structure of turbulence, namely that 
vortex structures dominate she e t-like structures resulting in an inhibition of the energy cascade 
( Avvalasomavaiula fc WarhaftJ ( 200fill l. 

The right hand side of Figure [5] shows the longitudinal kurtosis, a measure of the flow inter- 
mittency. Here the effects are milder in both directions, although a small increa se in the kurto¬ 


sis is noted in the expanding directions (from the expected value of about 5-6. iGvlfason et al. 
(I2004M . 


4 Particle acceleration statistics and discussion 

The non-uniform mean velocity field has a significant effect on the dynamics of tracers and 
inertial particles, both directly via the mean flow velocity and indirectly through the strained 
turbulent field. 
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4.1 Tracers accelerations in straining flow 

The full acceleration of tracer particles in our flow is given by the equation: 


Clpi 


Dup^ 

Dt 


Dup^ dUp, dUp^ 

Dt dxj dxj 


i= 1,2,3. 


( 11 ) 


The subscript p indicates that the full instantaneous flow velocity Up ., the fluctuating flow 
velocity Up ^, and the mean flow velocity Up. are taken at the location of the tracer. Here, we 
have assumed that the mean flow is time-independent. The mean of the tracer acceleration 
is Up. ; equal to the acceleration in a laminar flow of the same strain geometry. We are 
interested in the statistics of the tracer acceleration fluctuation, that is, when the pure mean 
flow contribution to the acceleration has been subtracted: 


^Pi 


Dup^ ^ ^Up^ 

Dt dx, 


i= 1,2,3. 


( 12 ) 


The first term on the right hand side refers to the material derivative of the fluctuating velocity 
held and represents the acceleration experienced by the huid particle advected by the huctuating 
how held, and the second term refers to the tracer acceleration induced by turbulent transport 
in the mean velocity held. 

In the variance of acceleration, the cross terms give rise to | ^Pi) (^® sum over i), 
representing the time derivative of the kinetic energy in the i-th component of velocity. In ad¬ 
dition, the latter term describes contributions of velocity variances, notably (25')^ ((^^pi)^| £md 

52 (('“^ 2)^1 for the hrst and second component, respectively. These terms are easily evaluated 
if the statistics of how velocity huctuations are available. 


4.1.1 Approximate tracer acceleration variances using rapid distortion theory 


When the straining is su fficiently rapi d, the nonlinear and viscous forces vanish from the Navier- 
Stokes equations (e.g. see lPond ( 200r)li i. and therefore their solution is particularly convenient in 
comparison to solving the full Navier-Stokes equations. Below, we re-derive the RDT predictions 
for the evolution of the huctuating velocity variances, as well as deriving the prediction of RDT 
on the evolution of the tracer acceleration variance. 

In RDT, each Fourier mode evolves independently. Let us consider a single mode of the 
huctuating velocity 

u{x,t) = 

The wave nu mber and the Fourier coefficients evolve according to the following set of equations 
( Pooel (l200nlP 



(13) 

(14) 


When the mean how geometry, U = {-2Sx,Sy,Sz), has been applied, equation (IT^ . results in 


Ki(t) = K°e 


(15) 
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where Si = -2S, S 2 = 
equation da, gives, 


S 3 = S, and l^p = + (^ 3)^6 



In particular, the long term prediction gives 


Similarly, for Uj, 


Ui{t) 


lO-2St 


U2{t) 


lO-St 


Usit) 




Taking the material derivative of the single mode flow fluctuating velocity, we have 


Dui 

~Dt 


-2Sui, 


Du 2 

~Dt 


-SU 2 , 


DU 3 

Dt 


-SU 3 . 


(16) 


Therefore the anisotropic contribution to {{DuijDt)"^) can be estimated as S]^{{uiY), where 
Si = -25', S 2 = S 3 = S are the strain rate at different directions. Together with the 
accele ration variances an = (l/3)(aiai)/((eo)^/j^)^^^ at the beginning of the straining 
(2002 ll l. the magnitude of tracer fluctuating acceleration variance can be approximated as 


normalized 
(Voth et al. 


((«Pi)^) “0 


(^) 


1/2 


+ 2Sf{{u,y) + S 


d{{u.Y) 

dt ’ 


*=1,2,3, 


(17) 


where the first term represents the isotropic contribution to acceleration variances with depen¬ 
dence on turbulence intensity, and all the terms on the right-hand side are Eulerian quantities 
of the straining flow. Here (cq) is the mean energy dissipation rate in homogeneous isotropic 
turbulence across all realizations, and (e) is the time-dependent mean energy dissipation rate 
in the straining turbulence across all realizations. 

Figure [5] shows the acceleration variances of passive tracers in the straining flow. The 
variance is higher in the compression direction than in the expanding directions due to the mean 
straining geometry, and the effects are seen immediately after the strain has been imposed. The 
approximations derived above, applying RDT, fit nicely to the simulation data as shown. The 
change in the isotropic dissipation rate due to straining is accounted for in the first term of 
da, and the rest of the terms involve the mean strain and the velocity fluctuations. The mean 
strain causes a rapid increase in the acceleration variance as the strain is applied, particularly 
at the higher rate of strain. The subsequent increase and the differences between the individual 
components are partially due to the evolution of the velocity variances for each component; 
namely that the compressed velocity components are emphasized (increasing energy content) 
whereas the expanding velocity components are either suppressed or maintained. 


4.2 Inertial particle accelerations in straining flow 


The differential equations for the particle position as a function of time, obtained by combining 
the mean flow components in equations m with the equations © and are second order 
linear ordinary differential equations with constant coefficients: 


dS'Xp. 1 dxp^ 

dt"^ Tp dt 



1 

— Ui, 


*=1,2,3. 


The roots Ai 2 = 


-1 ± + 'iSiTp 

2to 


of the characteristic equations of these equations prescribe 


the behavior of particle movements in absence of turbulence. In the compression direction, the 
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Figure 6: Normalized acceleration variances of passive tracers in the compressed and expanding 
directions vs. strain time at various strain rates. Solid and empty symbols represent the 
compressed and expanding directions, respectively. Solid (—) and dashed (—) lines indicate 
RDT long term predictions of the tracer acceleration variances from the straining together with 
normalized flow acceleration variances in HIT shown in (HZl). Dash-dot line (-.) indicates the 

term oq/ ((e)^/i^) for the S = 0.1 case. Diamond (O), circle (O), square (□), triangle (a) 
mark normalized tracer acceleration variances in strain rates S = 0.1,1,4 and 10. Left; tracer 
acceleration variances in different strain rates in R\q = 117. Right: tracer acceleration variances 
in S' = 4 in Rxq = 117 (squares (■,□): compressed and expanding directions, respectively) 
and 193 (asterisks (*) and pluses (+): compressed and expanding directions, respectively). 
An estimate of statistical error bar is shown and computed according to (fTIB with Xj being 
v in the j-th realization. 
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combination of strain rates S and the Stokes relaxation time Tp in our simulations give rise to 
the discriminant Di = 1 - SS'rp, which separates two possible movements: a overdamped-decay 
motion toward the stagnation plane a: = 0 or a underdamped-oscillation about the stagnation 
plane. In the expanding direction the discriminant Il 2,3 = 1 + 45'rp is always positive, and 
particles move away from the axis y = z = 0 exponentially in time when only mean flow is 
considered. The accelerations of particles follow a similar patten. 

When turbulent fluctuations are present in the strained flow, the statistical description of 
the dynamics of the inertial particles becomes much more complicated. Treating the fluctuating 
flow velocity as a source term, one can solve the equations formally using the Laplace transform. 
Specifically, the formal expression of particle accelerations are: 

In the compression direction: 


Opi(t) 


Ai - A; 


[(A^(Ai + 2S)e^^* - A?(A2 + 2S)e^^^)xpo + (A?e^^‘ - A^e^^*)up,o 


- r‘ui(a;p(r),T)(A?e^^(*-")-A^e^^(*-")) dr 

T„ Jo 


1 

Tp 


+ —Ui{Xp{t),t). 


(18) 


Hi < 0: 

api(<) = e”^ Upiocos(wt) +((1-4S'rp)upiO+ 8S'^rpa:po)sin(a;t) 


^ f ui(xp{T),T)e^ (cos{uj{t-T))+ — 3 :sin(w(t-T))'j 

ri Jo \ ZTpW / 


+ —Ui{Xp(t),t). 


(19) 


In the expressions. 


Ai, 2 = -;^(1 ± y/1 - SSTp), uj = -^y/8STp-l 

ZTp 


( 20 ) 


as well as the Stokes relaxation time Tp determine time scales for various contributions to inertial 
particle accelerations. Xpo and Up^o are position and velocity of the particle right before the 
straining starts (at t = 0“). In the expanding direction y: (the expression in the 2 direction is 
similar) 

“wW = -^[(A2(Ai-5)e^=*-A?(A2-^)e^^‘)ypo + (A?e^^*-A^e^^*)up,o 

Ai - A2 

+ A ru2(=rp(T),r)(A?e^^(*-")-A2e"=(‘-")) drl + lu2(=rp(t),t). (21) 

TV) -'0 Tri 


Similarly, ?/po and Vp^o are position and velocity of the particle at t = 0 , and 

Ai .2 =- 7^—(1 ± y/l + 45'rp). ( 22 ) 

ZTp 

The mean flow influences the variances when full acceleration is considered, and so is the case 
for other statistics involving full particle velocity or accelerations. Since the magnitude of the 
mean flow velocity depends on the location in the domain, the magnitudes in the variances of 
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Figure 7: Acceleration variances of inertial particles in the compressed (solid symbols) and 
expanding (empty symbols) directions vs. strain time for two types of particles Left: Tp = 0.015; 
Right and Tp = 0.05. Symbols refer to different strain rates, diamond (O) S = 0.1; circle (O) 
5=1; square (□) 5 = 4; triangle (a) 5 = 10. The data are from the data set R\o = 117. An 
estimate of statistical error bar is shown and computed according to cni) with Xj being particle 
acceleration variances in the j-th realization. 


acceleration components are characterized by the domain size in respective directions (through 
the initial positions XpOj VpO in acceleration expressions (IT^ . (I19|) and (EU) in addition to the 
rate of strain. 

In an attempt to minimize the influence of the mean flow, in addition to ensuring that 
our statistics are deduced from sufficient many independent samples, we condition our inertial 
particle analysis on particles that started in a thin layer parallel to and next to the x = 0 plane 
for the x-component statistics and a thin layer parallel to and next to the y = 0 plane for the 
y-component statistics. For the strain rates 5 = 0.1 and 5 = 1, we use layers with thickness of 
8 lattice units in the i?Ao = 117 simulations. For 5 = 4 and 5 = 10 flows, we use layers with 
thicknesses of 4 and 2 lattice units. With the selected thicknesses, we ensure the difference 
of the mean flow velocities within the layers do not exceed 55% of n™® in the compression 
direction, and do not exceed 77% of in the expanding direction. The numbers of particles 
available in the layers decreases with the thickness of the layers from about 61,500 in 5 = 1 
flows to 15,500 in 5 = 10 flows in the compression direction. In the expanding direction, the 
number of particles used for the analysis decreases from 240,000 in 5 = 1 to 60,700 in 5 = 10. 
We also apply the symmetry with respect to the planes a: = 0 and y = 0. 

Figure [7] shows the acceleration variances of inertial particles in the straining flow for two 
particle types Tp = 0.015 and Tp = 0.05, which correspond to St = 0.3 and 1 at the beginning of 
the straining. These two types of particles have positive discriminants in 5 = 1, correspond to 
L>i > 0 and Hi < 0 in 5 = 4, and have negative discriminants in 5 = 10. 
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Figure 8: Acceleration variances for tracers and inertial particles with Tp = 0.01,0.1. Lines with 
circles (O): Tp = 0.1; lines with triangles (a): Tp = 0.01. Solid lines with solid symbols: ((opj)^); 
dashed lines with empty symbols: ((opa)^)- Solid lines (—): tracers ((apj^); dashed lines (—): 
tracers {{ap 2 )^)- Left: A = 4, Right: S = 10. The data are from the data set R\q = 117. An 
estimate of statistical error bar is shown and computed according to (fTUl) with Xj being particle 
acceleration variances in the j-th realization. 


4.2.1 Initial transition period of acceleration variances 

The acceleration variances of the inertial particles show a transition period at the beginning of 
the straining that is not seen in the tracer accelerations. Equations m and (USD indicate that 
the initial position Xpo and velocity Vp-^o of an inertial particle affects its acceleration. The time 
scale of this influence depends on the exponents of the exponential terms in the acceleration 
expressions. In the compression direction, if Hi < 0 the decaying rate is 2rp. If Di > 0, both 
Ai and A 2 are negative so the decaying rate is min(|^, |^) = p^. This explains the different 
lengths of the initial transition period for particles with different Stokes numbers in flows with a 
hxed strain rate. In the left figure of Figure HI we show acceleration variances of particles with 
Tp = O.OI and O.I in S' = 4 to elucidate the time scale for initial transition. The ratio between 
the decaying rates of these two types of particles is about 1.75. 

Although (( 1 x 2 )^) decreases with straining, ( 0 ^ 2 ) increases at the beginning and a maximum 
occurs. This is because A 2 is positive in the expanding direction. Similar to the case discussed in 
the cc-direction, the exponents determine the time scale of the behavior of acceleration variances. 
From (EU one can estimate the time that the maximum takes place by ignoring the flow 
fluctuation, and obtain the ratio between the times of maximum acceleration variances of the 
two particles presented in the figure to be 1.52. It appears that our simulation data confirms 
such an estimation. For S = 10, our simulations were too short to display the post-transition 
period. 

4.2.2 Magnitude of the acceleration variances 

Figure [7] shows that the acceleration variances of inertial particles increases with the strain rate. 
This is mainly due to the mean flow contributing to the acceleration in terms of the factors 2S 
and S in the expressions (fTSl) . (fTOl) and (ITT!) and, in a smaller part, through the exponents that 
regulate the influence of the initial conditions. 
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At the onset of straining, the acceleration variances in the expanding direction is higher 
than that in the compression direction. This could be explained through the exponents of the 
kernels. In the compression direction, the coefficients Ai,A 2 and in the exponents of the 


kernels in (fl^ and (HU are all negative and indicate decay with time. But for the expanding 
direction, A 2 = jg positive, which leads to the increase of magnitude. For sufficiently 

long time, the influence of the initial condition is diminished. The magnitude of the flow velocity 
fluctuation (in the last two terms of the acceleration expressions) in the expanding direction 
remains roughly constant, while in the compression direction the flow velocity fluctuation grows 
exponentially. As a result the acceleration variances in the compression direction overtakes that 
in the expanding direction. 

In contrast to the isotropic homogeneous turbulence ( Avvalasomavaiula et al. ( 2006ll : lBec et al\ 
( 2006ll b the acceleration variances of inertial particles do not necessary have lower magnitudes 
compared with that of tracers in strained flow. The inertial particle with Tp = O.OI in the 
right plot of Figure [5] depicts such a situation. Recall that (flTll provides an approximation of 
tracer acceleration variances, and in higher strain rates the terms 2 Sf{{ui)^) + Sid{{ui)^)Idt 
are the main contributor of the variances. From Figure [3] we know that the variance of flow 
fluctuating velocity grows exponentially in time, so its time derivative can be estimated as a 
constant multiple ciS (ci is a constant) of the variance itself. Hence the main terms in the 
tracer acceleration variance increases in the order of Sf{{ui)^) with the strain rate. For the 
inertial particles, however, their long term variances depends on the term {{ui)'^)lTp. When Tp 
is small enough, the factor I/r^ is larger than the magnitude of Sf, and the inertial particle 
acceleration variances surpass the tracer acceleration variances. 


4.3 Probability density functions of particle accelerations 


Figures [5] to [TT] show the particle acceleration pdfs at various rates of strain {S = 1,4 and 10) and 
at two Reynolds numbers {Rxo = 117 and R\o = 193) for tracers and inertial particle {Tp = 0.015 
and Tp = 0.05), respectively. The effect of increasing the rate of strain is most notably seen by 
the narrowed pdf-tails for the tracer and inertial particle accelerations. 

The top left plot of Figure IHl demonstrates the narrowing of the tracer acceleration pdfs 
in straining flow. The top right plot shows that the narrowing effect is milder at the higher 
Reynolds number due to the faster time response of the smaller scales at the higher Reynolds 
number. The increase in the magnitude of acceleration variances due to the mean straining 
appears to be the primary reason for the tail-narrowing. The bottom left plot of Figure [3] 
indicates the response of acceleration in the compressed and expanding directions. Although 
the flow field is very different component wise, the acceleration pdfs, resulting from tracers 
following trajectories of small scale structures in the fluid are more or less identical. Addition¬ 
ally, the evolution of the acceleration flatness ((op^ - Op.) )/((api - OpJ (a measure of the 
intermittency of the acceleration) is shown in the bottom right figure for the lower Reynolds 
number R\q = 117, emphasizing the narrowing effect of straining and the difference between 
components. 

In terms of inertial particles, as for isotropic homogeneous turbulence, the narrowing of the 
pdf tails follows an increase in the Stokes number. Figure flOl illustrat es the narrowing effects. 
Such behavior ha s been demonstrated in a number of previous studies! Avvalasomavaiula et al\ 
(l2006^ : lBec fJ^(I 2006IB . owing to the heavier particles passing through, or selectively filtering, 
the most rapid motions in the flow field. However, here the situation is more complex, since 
the acceleration variance is not necessarily smaller for inertial particles, due to the strong effect 
of the mean flow as discussed above. Figure [TT] shows that the higher Reynolds number has 
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Figure 9: Normalized pdfs of tracer accelerations and tracer acceleration flatness. Top left: 
Pdf(apj) in HIT (circle •), straining flow 5” = 4 (square ■), and 5” = 10 (triangle A). Top 
right: Pdf(api) in HIT (circle •: R\o = 117, o : = 193) and in straining flow S' = 10 

(trianlge A : i?Ao = 117, A : _Rao = 193). Bottom left: Pdf(apj) and Pdf(ap 2 ) in Rxo = 117 
flow. Circles (•, o): HIT, squares (■,□): S = 4, triangles (A, A): S = 10. Solid sym¬ 
bols: i = 1 component, empty symbols: i = 2 component. Bottom right: Tracer acceleration 
flatness,((flp. -a^Y)l{{a,pi at R\q = 117. Solid line (—): HIT, squares (■,□): S = 4, 

triangles (A, A): S = 10. Solid symbols: i = \ component, empty symbols: i = 2 component. 
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Figure 10: Normalized probability distribution function of particle acceleration components ap. 
for different strain rates, for particles originating in thin slices parallel to the a; = 0 plane, for the 
i = 1 component and to the y = 0 plane for the i = 2 component. Left: particles with St^ = 0.3 
in R\a = 117. Right: particles with St^ = 1 in R\q = 117. Top: The i = 1 component. Bottom: 
The i = 2 component. Circles (•): HIT; squares (■): S' = 4; triangles (a): S = 10. 


a expanding effect on lighter inertial particle acceleration pdf tails. However, for the heavier 
particles {St = 1) in higher strain S = 10, the effect is not as evident in the expanding direction. 
We note that the Stokes number of a particle increases slightly during the straining due to 
a decrease in the Kolmogorov timescale (less than 10% for the highest rate of strain). This 
effect contributes to the narrowing of the acceleration pdf tails, but to a lesser effect than the 
increased variances. 

It is interesting to consider our results in a wider context of flows with non-zero mean 
co mponents, for exampl e by co mparing with the dynamics in shear flow. Lavezzo et al. ( 20I0ll 


and iGerashchenco et al\ (12001^ considered tracers and inertial particles in a non-uniform shear, 
namely a turbulent boundary layer. An increased rate of shear, closer to the wall, resulted in 
an increased acceleration variance in the lateral component but a milder effect in the transverse 
component, which appears to be consistent with the predictions of equation m for their 
geometry. As a result of the increased acceleration variances, attenuation was found in the tails 
of the inertial particle acceleration pdfs as observed in this work. 
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Figure 11: Normalized probability distribution function of particle acceleration components ap. 
for two Reynolds numbers, for particles originating in thin slices parallel to the a; = 0 plane, 
for the 1=1 component and to the y = 0 plane for the i = 2 component. Left: particles with 
Stf) = 0.3 in R\o = 117 and Sto = 0.34 in Rxq = 193. Right: particles with Sto = 1 in R\o = 117 
and Sto = 1.12 in R\q = 193. Top: The i= \ component. Bottom: The i = 2 component. Circles 
(•, o): HIT; triangles (A, A): S = 10. Solid symbols: i?Ao = 117; empty symbols: i?Ao = 193. 


20 






5 Conclusions 


Our results show a strong effect of the mean flow straining on the Lagrangian acceleration 
statistics, both for the passive tracers and the inertial particles. The effect of straining is 
primarily felt in acceleration variances and pdfs when the rate of strain is sufficiently high such 
that the strain timescale is of a comparable magnitude to the Kolmogorov timescale. For high 
rates of strain the magnitude of the acceleration variances are increased significantly and the 
tails of the normalized acceleration pdfs for tracers and inertial particles narrow. The former 
effect is well explained by observing the predicted behavior of the acceleration variance by 
rapid distortion theory. RDT provides us a relation between the flow velocity variances and its 
acceleration variances and illustrates the dependence of the acceleration variance on the rate of 
strain. 

However, the effect is complex, partly due to the connection, or lack thereof, between particle 
acceleration component in one direction and the fluid flow in the same direction; a particle 
trajectory around a strong vortex will result in large acceleration values in the directions normal 
to the axis of the vortex. But there is also a direct contribution from the mean straining and 
fluctuating velocity in the acceleration, resulting for example in an increased variance value of 
acceleration in both compressing and expanding directions. 

For tracers, the narrowing of the normalized acceleration pdfs stems therefore in a complex 
manner from both of these effects. The same effect is also felt by the lighter inertial particles, 
given that their inertia is sufficiently small to sample the small scale motions. Because of their 
small inertia the interplay with the mean flow enable their acceleration variances to rise beyond 
that of tracers. 

When the inertia is further increased the ballistic particle motion in the rapidly accelerating 
mean flow becomes increasingly important leading to Gaussian pdf tails. Here, the particles 
are swept through the fluid, and the slower, large scales of motions are more likely to influence 
the particles. This could also be seen by the lower magnitude of acceleration variances of the 
heavier particles compared with the lighter inertial particles. We derive the formal expressions 
for inertial particle acceleration, and these expressions reveal the complex interplay between 
flow straining and particle inertia. 

Our findings emphasize the importance of the presence of strong mean motions and imposed 
small scale anisotropy in particle laden flows. It is our opinion that the results have relevance 
to the understanding and modelling of a range of practical deforming or straining flows where 
inertial particles are important aspects of the process. In particular, we believe that our findings 
may help in the development of sub-grid Lagrangian models for particles in the proximity of 
straining regime near stagnation points. 
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